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Abstract 


Application  of  Total  Variation  Diminishing  (TVD)  schemes  to  turbulent  flows 
is  considered.  The  mathematical  and  physical  basis  of  TVD  schemes  is  discussed. 

TVD  methodology  is  extended  to  the  solution  of  turbulent  flow  problems.  A 
first-order  time  accurate,  second-order  space  accurate  algorithm  is  used  to  compute 
solutions  to  the  problems  of  shock-boundary-layer  interaction,  turbine  rotor  cascade 
flow,  and  unsteady,  shock-induced  heat  transfer  using  the  TVD  algorithm.  This  al¬ 
gorithm  provides  the  capability  to  accurately  predict  separation,  reattachment,  and 
pressure  and  skin  friction  profiles  for  shock-boundary-layer  interaction.  Improved 
accuracy  is  demonstrated  in  computing  surface  pressures  for  a  turbine  rotor  cas¬ 
cade.  Heat  transfer  for  the  cascade  is  predicted  with  fair  accuracy,  showing  aU  the 
significant  features  of  the  experimental  Stanton  number  profile.  Fairly  accurate  com- 
pcirison  with  theory  and  experiment  is  evident  for  the  unsteady,  shock-induced  heat 
transfer  problem,  with  the  exception  being  failure  to  correctly  predict  transition  for 
this  case. 
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IMPROVED  MODELING  OF  UNSTEADY  HEAT  TRANSFER 

(The  First  Step) 


/.  Mathematical  and  Physical  Basis  of  TVD  Schemes 

1.1  The  Genesis  of  TVD 

Total  Variation  Diminishing  (TVD)  schemes,  originally  refered  to  as  Total  Vari¬ 
ation  Nonincreasing  (TVNI),  first  appeared  in  1983  with  the  pubhcation  of  Harten’s 
High  Resolution  Schemes  for  Hyperbolic  Conservation  Laws  [23].  In  general,  TVD 
schemes  are  arrived  at  by  applying  a  first-order  accurate  numerical  method  to  an 
appropriately  modified  flux  function  thus  yielding  a  method  that  is  second-order  ac¬ 
curate  except  near  points  of  extrema  of  the  solution.  The  genesis  of  the  TVD  class 
of  finite-difference  schemes  can  be  traced  to  1976  when  Harten,  Hyman,  and  Lax  au¬ 
thored  On  Finite- Difference  Approximations  and  Entropy  Conditions  for  Shocks  [24]. 
This  work  first  addressed  the  question  of  whether  finite-difference  approximations 
to  the  solution  of  hyperbolic  conserration  laws  converge  to  the  physically  relevant 
solution.  This  is  of  interest  because  weak  solutions  to  such  conservation  laws  are  not 
uniquely  determined  by  initial  values,  but  require  an  entropy  condition  be  met  to 
converge  to  the  particular  physical  solution  [24]. 

In  the  mid  1970’s,  Harten  was  also  working  on  his  artificial  compression  method 
(ACM)  [20]  to  modify  standard  finite-difference  schemes  in  an  effort  to  prevent  the 
smearing  of  contact  surfaces  and  improve  shock  resolution  [21,  22].  Prior  to  this 
effort,  Harten  states  that  the  standard  finite-difference  schemes  in  use  typically 
smeared  shocks  over  3-5  cells  while  the  width  of  the  contact  surface  behaved  as 
^i/(H+i)  n  is  the  total  number  of  time  steps  taken  and  R  is  the  scheme’s  order 

of  accuracy.  Harten’s  ACM  also  addressed  the  fact  that  schemes  of  order  greater 
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than  one  produced  overshoots  and  undershoots  around  the  discontinuity  [21]  and 
forced  the  approximated  solution  to  be  nonphysical  [24].  Harten’s  ACM  modifica¬ 
tions  to  existing  schemes  provided  the  foundation  for  the  new  class  of  TVD  schemes 
presented  in  his  1983  paper. 

The  rigorous  mathematical  foundation  of  TVD  schemes  is  mainly  confined  to 
scalar  hnear  and  nonlinear  conservation  laws  and  is  painstakingly  outhned  in  ref¬ 
erences  [24]  and  [23].  Computational  fluid  dynamicists  are  interested  in  applying 
TVD  schemes  to  systems  of  nonhnecir  hyperbohc  conservation  laws,  such  as  the  Eu¬ 
ler  equations  of  gasdynamics.  Therefore,  Harten  details  the  application  of  TVD 
methodology  to  1-D  systems  using  Roe’s  approximate  Riemann  solver  and  provides 
an  example  of  its  extension  to  2-D  using  Strang’s  dimensional  splitting  [23].  The  orig¬ 
inal  Harten  scheme  was  a  second-order  accurate  exphcit  method  but  was  extended 
to  a  second-order  accurate  imphcit  method  by  Yee  and  Harten  [40]. 

The  high-resolution  TVD  approach  soon  gathered  favor;  explicit  and  implicit 
variations  were  then  applied  to  the  Euler  equations  in  general  geometries  by  Yee  and 
Kutler  [41]  and  by  Yee  and  Harten  [43].  Later,  Wang  and  Widhopf  further  extended 
Harten’s  TVD  methodology  to  a  finite- volume  scheme  for  the  Euler  equations  [44]. 
TVD  algorithms  have  continued  to  develop  over  the  past  decade.  Harten’s  origi¬ 
nal  scheme  was  of  the  upwind  Vciriety,  meaning  that  the  modifications  to  the  flux 
function  are  apphed  based  on  the  direction  of  wave  propagation,  or  characteristic 
direction.  Symmetric  algorithms  have  since  come  into  use  where  the  modifications 
are  applied  without  regard  to  the  characteristic  directions.  Methods  are  also  avail¬ 
able  for  partial  differential  equations  with  source  terms  and  stiff  source  terms.  Yee’s 
1989  publication, A  Class  of  High- Resolution  Explicit  and  Implicit  Shock-Capturing 
Methods  [42],  provides  detailed  information  on  numerous  versions  of  TVD  algorithms 
and  examples  of  their  application  to  numerous  problems. 
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1.2  Hyperbolic  Conservation  Laws  and  TVD  Methodology 

The  present  section  provides  a  description  of  the  hyperbolic  conservation  laws 
for  which  TVD  schemes  provide  solutions.  The  requirements  for  uniqueness  of  a 
solution  to  the  initial  value  problem  are  given  along  with  the  necessary  conditions 
to  guarantee  convergence  of  a  finite  difference  approximation  to  this  solution.  A 
summary  is  provided  of  the  methodology  behind  the  construction  of  Harten’s  original 
second-order  accurate  TVD  scheme. 

1.2.1  Finite- Difference  Schemes  and  Oleinik’s  Entropy  Condition.  The  present 
analysis  is  concerned  with  weak  solutions  of  the  initial  value  problem 

ut  +  f{u)x  =  0 

—  OO  <  X  <  CO  (1.1) 

u(x,  0)  =  <p(x) 

where  u{x,t)  is  a  column  vector  of  m  unknowns,  f{u)  is  the  flux  vector  of  m 
components,  and  (f){x)  is  the  initial  data.  Eq  1.1  is  hyperbohc  if  all  eigenvalues 
a^(u),...,a”'(u)  of  the  Jacobian  matrix 

A(,u)  =  U  (1.2) 

are  real  and  the  set  of  right  eigenvectors  /?*(«),...,  ft™  (a)  is  complete  [23]  over  the 
domain. 

Following  Harten  [23],  consider  systems  of  conservation  laws,  Eq  1.1,  possess¬ 
ing  an  entropy  function  U (u)  defined  such  that 

>  0 

t/u/„  - 

where  F  is  a  function  known  as  the  entropy  flux  [23]. 


(1.3) 


3 


The  class  of  all  weak  solutions  to  Eq  1.1  is  too  laxge  in  that  the  initial  value 
problem  is  not  unique  [24].  An  additional  constraining  relation  is  needed  if  the 
scheme  is  to  choose  the  physicEiUy  relevant  solution.  This  additional  constraint  is 
known  as  Oleinik’s  entropy  condition  and  can  be  expressed  as  [24] 

Uiu)t  +  F{u),<0  (1.4) 

Let  us  now  consider  numerical  solutions  to  Eq  1.1  obtained  using  a  {2k  +  1) 
point  explicit  scheme  in  conservation  form  [24].  A  scheme  is  in  conservation  form  if 
it  can  be  expressed  as 


»r' = ^  (1-6) 

where 

/J+:/2  =  / 

and  A  =  At  I  Ax.  In  Eqs  1.5  and  1.6,  /  is  the  “numerical”  ,  or  mesh,  flux  function 
consistent  with  f{u)  in  that  f{u, . . .  ,u)  =  f{u).  The  solution  u  is  approximated  on 
the  mesh  by  =  v{jAx,nAt).  The  numerical  scheme  given  by  Eq  1.5  is  consistent 
with  the  entropy  condition,  Eq  1.4,  if 

<  u;  -  A  (i.r) 

where  Uf  =  U  [vj),  =  F  ,  md  F  is  the  numerical  entropy 

flux  consistent  with  F{u)  such  that  F{u, . . .  ,u)  =  F{u) 

The  question  of  convergence  of  the  finite  difference  scheme,  Eq  1.5,  to  the  ap¬ 
propriate  weak  solution  of  Eq  1.1  must  now  be  addressed.  The  scheme  under  consid¬ 
eration  is  nonlinear,  so  stabiUty  of  a  consistent  scheme  does  not  imply  convergence. 
Harten  [23]  outlines  three  conditions  which,  when  satisfied,  ensure  convergence. 
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(1)  The  total  variation  {TV)  of  the  finite  difference  scheme  is  uniformly 
bounded,  where 

OO 

TV(v)=  £  (1-8) 

j--oo 

(2)  The  scheme  is  consistent,  as  Ax  — >  0  ,  with  Oleinik’s  entropy 
condition  for  aU  entropy  functions  of  Eq  1.1. 

(3)  Oleinik’s  entropy  condition  impUes  a  unique  solution  of  the 
initial  value  problem  for  Eq  1.1. 

The  reader  is  referred  to  the  references  given  by  Harten  [23]  for  the  arguments 
that  imply  convergence  given  satisfaction  of  the  above  criteria.  For  the  present 
work,  the  validity  of  these  criteria  will  be  assumed  and  the  effort  concentrated  on 
demonstrating  the  development  of  a  scheme  that  satisfies  criteria  (1)  and  (2)  when 
given  the  third  criterion. 

J.S.S  Development  of  Marten’s  Second- Order  Scalar  TVD  Scheme.  Marten’s 
second-order  accurate  TVD  scheme  is  the  product  of  a  nonoscillatory,  first-order 
accurate  scheme  applied  to  an  appropriately  modified  flux  function  [23].  This  section 
describes  the  properties  of  the  first-order  scheme  and  outhnes  the  procedure  used  by 
Harten  to  mrive  at  the  appropriate  modified  flux. 

Consider  the  initial  value  problem  for  a  scalar  conservation  law: 

Ut  +  f{n)x  =  +  a{u)u^  =  0 

a{u)  =  /u  -OO  <  2  <  oo  (1-9) 

u{x,  0)  =  (j){x) 

where  (f){x)  is  of  bounded  total  variation.  Rigorous  analysis  is  restricted  to  the  sceilar 
case  because  TVD  schemes  are  not  defined  for  systems  of  of  non-hnear  conservation 
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laws  where  the  spatial  total  variation  of  the  solution  may  increase  due  to  wave 
interaction  [23]. 

A  weak  solution  of  Eq  1.9  has  a  monotonicity  property  [23],  as  a  function  of 
time,  defined  as: 

(1)  No  new  local  extrema  in  x  may  be  created. 

(2)  A  local  minimum  is  nondecreasing  and  a  local  maximum 
is  nonincreasing. 

The  monotonicity  property  impHes  that  the  total  variation  in  x  is  nonincreasing  in 
time,  TV{u{t2))  <TV{u{t^)). 

An  exphcit,  {2k  +  1),  point  finite-difference  scheme  in  conservation  form,  as 
given  by  Eq  1.5  and  applied  to  Eq  1.9,  can  be  written  as 

“r  =  . »?+*)  _  (1,10) 

or  in  operator  notation  as 

=  (1.11) 

The  scheme  given  by  Eq  1.10  is  TVD  if,  for  all  v  of  bounded  total  variation 

TV{L-v)  <TV{v)  (1.12) 

where  the  total  variation  is  defined  by  Eq  1.8.  Eq  1.11  represents  a  monotonicity 

preserving  scheme  if  the  operator  L  is  monotonicity  preserving.  That  is,  if  u  is  a 
monotonic  mesh  function  so  is  L  •  v.  The  scheme  is  monotone  if  is  a  monotonic 
nondecreasing  function  of  each  of  its  2fc  -f  1  arguments  [24]: 

>0  (1.13) 
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for  all  i  such  that  —k  <  i  <  k. 

An  example  of  a  scheme  that  is  not  monotone  is  the  second-order  accurate 
Lax-Wendroff  scheme  with 


/i+i  =  5  [/  (•’,“)  +  f  (“"+.)  - 

where  ^j+LV  —  Uj+i  —  Vj  .  Therefore  the  discrete  equation  is 

IM/ (”"+.) --f  2*-? +”?->)]  (1.15, 

Taking  the  derivative  of  H  with  respect  to  the  argument  ■^i+i  yields 

=  5  -  ■')  (116) 

where  v  =  a\  .  Only  the  case  0  <  <  1  need  be  examined  since  the  Law-Wendroff 

scheme  is  unstable  for  z/  >  1,  and  Lax-Wendroff  provides  the  exact  solution  for  z/  =  1. 
Clearly,  the  Lax-Wendroff  scheme  is  not  monotone  for  any  i/  <  1  .  Additionally, 
the  numerical  results  of  reference  [24]  show  that  the  Lax-Wendroff  scheme  is  not 
monotonicity  preserving. 

The  first-order  accurate  Roe  scheme  provides  an  example  of  monotone  behav¬ 
ior.  The  numerical  flux  for  the  Roe  scheme  is 

=  i  [/  (»?)  +  /  (1.17) 

giving  the  discrete  equation 

(,7,, -2.7 + 

=  7f(»7_„»7,^7^,) 
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Taking  derivatives  of  H  with  respect  to  each  of  its  arguments  gives 

^ 

H„.  =  1-u  (1-19) 

H.,..  =  0 

Thus,  if  is  a  monotonic,  non-decreasing  function  of  each  of  its  arguments  showing 
that  the  Roe  scheme  is  indeed  monotonic. 

Let  Sm,  Stvd,  and  Smp  denote  monotone,  TVD,  and  monotonicity  preserving 
schemes  respectively.  Theorem  2.1  of  reference  [23]  provides  the  hierarchy  of  these 
properties: 

Sm.  C  Stvd  C  Smp  (1-20) 

Thus,  the  Roe  scheme  is  also  TVD  and  monotonicity  preserving. 

A  scheme  in  the  conservation  form  of  Eq  1.10  that  is  monotone  with  converg¬ 
ing  boundedly  almost  everywhere  to  some  function  u[x,t)  has  two  further  desirable 
properties.  The  theorem  of  Lax  and  Wendroff  as  given  by  reference  [24]  states  that 
if  the  scheme  is  in  conservation  form  with  v{x,t)  converging  almost  everwhere  to 
u{x,t),  then  u{x,t)  is  a  weak  solution  of  Eq  1.9.  The  theorem  of  Harten,  Hyman, 
and  Lax  [24]  states  that  if  the  scheme  is  monotone  in  addition  to  meeting  the  crite¬ 
ria  of  the  Lax- Wendroff  theorem,  then  Oleinik’s  entropy  condition  is  satisfied  for  aU 
discontinuities  of  u.  Thus  a  monotone  scheme  satisfies  the  convergence  criteria  for  a 
unique  solution  of  the  initial  value  problem  as  stated  in  the  Section  1.2.1. 

Attention  is  now  focused  on  how  the  properties  of  a  monotone  scheme  are  help¬ 
ful  in  constructing  Harten’s  second-order  TVD  scheme.  Harten  states  that  monotone 
schemes  provide  second-order  accurate  solutions  to  the  modified  Eq  [23] 

ut  + f{u)^  =  At\/3{u,X)ux]^  (1-21) 
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(1.22) 


/?(u,A)  = 


2A2 


l=-k 


P[u,X)  >  0 


/3{u,  A)  7^  0 

where  j3  is  a.  numerical  dissipation  term.  Since  j3[u,X)  ^  0,  monotone  schemes  are 
only  first-order  accurate  approximations  to  the  initial  value  problem  of  Eq  1.9. 

Suppose  the  scheme  given  by  Eq  1.10  is  a  monotone  scheme  and  thus  provides 
a  second-order  accurate  numerical  approximation  to  the  modified  equation,  Eq  1.21, 
rewritten  as 

+  =  0  (1-23) 

where  g  =  AxP{u,X)u^.  Applying  this  scheme  to  the  following  equation 


Ut  +  if  +  {^/^)9)x  — 


(1.24) 


yields  a  second-order  accurate  approximation  to  its  modified  equation.  Since  g  =  (9 [As] 
the  modified  equation  satisfies  [23] 

„,  +  /.  =  0  [(A^;)^]  (1.25) 


Thus,  application  of  a  first-order  scheme  to  a  scalar  conservation  law  with  an 
appropriately  modified  fiux  function  yields  a  second-order  accurate  approximation 
to  the  original  equation  Ut  +  fx  =  0.  Note  that  in  order  to  apply  the  scheme  to  the 
modified  flux  function,  g  must  be  a  differentiable  function  of  u.  Harten  achieves  this 
by  smoothing  the  point  values  of  g  [23].  This  smoothing  enlarges  the  support  of 
the  scheme  such  that  his  first-order  scheme,  using  a  three-point  stencil,  becomes  a 
second-order  scheme  using  a  five-point  stencil.  The  reader  is  referred  to  reference  [23] 
for  the  details  of  how  the  three-point,  first-order  scheme  is  constructed  so  as  to  ensure 
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its  TVD  property. 

Let  us  now  turn  our  attention  to  the  specific  scalar  scheme  developed  by 
Haxten.  Consider  a  three-point,  finite-difference  scheme  in  conservation  form  with 
the  following  numerical  flux  function 

/ {vj,vj+i)  =  2  [-^  ^  (^i+i)  "  {^/^)Q  (-^S+i/z)  ^j+i/2^^]  (1-26) 

where  Aj^i/2'^  —  '^j+i  ~ 

a  =  [/K-+i) -/(^i)]Mi+i/2^  (Ai+i/2i^  7^  0) 

=  aivj)  (^i+1/22^  =  0) 

Q  is  a  function  known  as  the  coefficient  of  numerical  viscosity.  Numerical  viscosity 
is  the  mechanism  that  allows  a  discontinuity  to  be  captured  as  part  of  the  numerical 
solution  [21].  This  is  in  contrast  to  shock  fitting,  where  the  discontinuity  is  considered 
as  an  internal  boundary. 

Lemma  3.1  of  reference  [23]  states  that  the  above  scheme  is  TVD  under  the 
Courant-Friedrichs-Lewy  (CFL)  condition 


A  max 
j 


< 


(1.28) 


given 

lx|  <  Q{x)  <  1  (1.29) 

for  0  <  la:|  <  /i  <  1. 

The  first-order  accurate  three-point  scheme  given  by  Eq  1.26  is  converted  to 
a  second-order  accurate  scheme  by  applying  the  three-point  scheme  to  modified  flux 
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values  [23]: 

fj^  =  /  {vj)  +  9j  ~  9 

^i+i/2  =  ^i+i/2  +  7i+i/2  7i+i/2  =  (fl'i+i  -  9j)  l^j+il^v 

where  u  =  Aa.  The  modified  numerical  flux  is 

=  I  [/f  +  -  (lAW 

=  |[/(^j)  +/(^i+l)] 

T  (l/(^'^))  [^i  9j+i  ~~  Q  (^i+i/2  T  ^i+i/z”^] 


(1.30) 


(1.31) 


Lemma  3.2  of  reference  [23]  provides  that  Eq  1.31  represents  the  numerical  flux 
of  a  second-order  scheme  so  long  as  Q(a;)  is  Lipschitz  continuous  and  gj  satisfies 


9j  +  5^1+ 1 

7j+l/2^i+l/2'y 


Q  (^^i+l/2)  -  (^i+1/2)  ]  '^j+1/2^  +  O  [2] 

9Hi-9^  =  0{A^) 


(1.32) 


Harten  [23]  constructs  g  in  the  following  manner  so  as  to  satisfy  Eq  1.32: 


9j 


•Sj+1/2  max  0,  min  (  gj+1/2  ,9j-i/2  ■  ■5i+i/2)] 


=  ^i+1/2  mm 

=  0 


9j+l/2 


9j-l/2\ 


\  (1-33) 

(ffi+1/2  •  ^i-1/2  >  Oj 
(^7  +  1/2  •  9i-\l2  <  0) 


where  ^ 

57+1/2  =  i  [Q  (^7+1/2)  -  (^7+1/2) 

"S7+I/2  =  -55^  {9j+l/2) 

Finally,  Lemma  3.4  of  reference  [23]  provides  that  a  conservative  finite  differ¬ 
ence  scheme,  with  the  numerical  flux  given  by  Eq  1.26,  is  TVD  under  the  restriction 
of  Eq  1.28  so  long  as  Q{x)  satisfies  Eq  1.29.  Thus  a  second-order  accurate  (except 
near  points  of  extrema  where  Sj^if2  is  discontinous),  five-point  scheme  has  been  con- 


^7+1/2^ 


(1.34) 


11 


structed  for  the  solution  of  Eq  1.9.  The  scheme  provides  high  resolution  capturing 
of  discontinuities  and  converges  to  a  physically  relevant  solution. 

1.2.3  Extension  to  Systems  of  Conservation  Laws.  We  now  concern  ourselves 
with  extending  the  scalar  scheme  developed  in  Section  1.2.2  to  systems  of  conserva¬ 
tion  laws.  Currently,  TVD  schemes  are  only  defined  for  scalar  hyperbohc  conserva¬ 
tion  laws  or  constant  coefficient  hyperbolic  systems.  This  is  due  to  the  fact  that  the 
spatial  total  variation  of  the  solution  to  a  system  of  nonhnear  conservation  laws  is 
not  necessarily  a  monotonically  decreasing  function  of  time  [43].  Wave  interactions 
may  cause  the  total  variation  to  increase.  Harten  extends  the  technique  using  a  gen- 
erahzed  version  of  Roe’s  approximate  Riemann  solver  [23].  The  idea  is  to  apply  the 
scheme  in  a  scalar  fashion  to  each  of  the  systems  hnearized  characteristic  variables. 

After  Harten  [23],  let 


S{u)={R'{u),...,R'"(u))  (1.35) 

be  a  matrix  whose  columns  are  the  right  eigenvectors  of  the  Jacobian  matrix  A{u) 
in  Eq  1.1.  It  follows  that 

S-^AS  =  A  (1.36) 

where  A  is  the  diagonal  matrix  of  eigenvalues  such  that  A,j  =  a'{u)Sij.  Therefore 

-f  =  0  (1.37) 

or 

+  A5-'ux  =  0  (1.38) 

where  the  characteristic  variables  w  axe  defined  such  that 

w  =  S~^u  (1.39) 
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Eq  1.38  becomes 


Wt  +  A.Wx  =  0 


(1.40) 


which  can  be  decoupled  into  m  sccJar  chciracteristic  equations  with  1  <  A:  <  m 

+  a^wl  =  0  (1.41) 

The  most  beneficial  use  of  the  chciracteristic  variables  comes  to  fight  by  rec¬ 
ognizing  that  they  can  be  viewed  as  the  components  of  u  in  the  coordinate  system 
such  that  [23] 

m 

(1.42) 

fc=i 

Harten  uses  this  fact  to  extend  his  scalar  scheme  to  general  nonlinear  systems  of 
hyperbolic  conservation  laws. 

Let  ct^+i/2  be  the  component  of  Aj4.i/2'y  =  Vj+i  —  'Vj  in  the  {R'^}  coordinate 
system  such  that 

m 

Aj4-i/2n  -  Q^+i/2^t+l/2  (1-43) 

Jt=l 

The  scheme  given  by  Eqs  1.30-1.34  is  extended  to  general  systems  as 


(1.44) 

/i+i/2  ^ 

=  l[/(^i) +  /(^i+l)] 

(1.45) 

where 

=  '^“J+1/2  and 

3)  =  ^J+1/2  max  [0, min  (I5J+1/2I  > 5J-1/2  ’  ^5+1/2)] 

(1.46) 
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with 

5j+l/2  =  2  (^i+1/2)  “  ('^i+1/2)  ]  “i+1/2 

^i+1/2  =  ^9n{9j+i/2) 

7-+1/2  =  (^i+l  -  ^?)  /«i+l/2  (“i+l/2  7^0) 

=  0  (“i+1/2  =  0) 


(1.47) 


1.2.4  Entropy  Enforcement.  As  a  final  comment  on  the  initial  development 
of  Harten’s  second-order  TVD  scheme,  we  turn  now  to  the  question  of  physically 
relevant  solutions  for  systems  of  equations.  As  mentioned  in  the  previous  section, 
the  total  variation  may  not  be  a  monotonic  decreasing  function  of  time  due  to  wave 
interactions.  In  addition,  Oleinik’s  entropy  inequahty  ensures  physically  relevant,  or 
admissable,  solutions  only  in  the  Hmit  as  Aa;  ->  0.  In  reahty  we  axe  concerned  with 
obtaining  admissable  solutions  on  a  relatively  coarse  mesh. 

In  order  to  arrive  at  a  proper  criterion,  Harten  examines  the  Riemann  initial 
value  problem  [23]  for  Eq  1.1: 

u(a:,0)  =  f(x)  =  Ul  x<0 

=  Ur  X  >  0 

with  Ur  and  ur  satisfying  the  Rankine-Hugoniot  relations  with  wave  speed  s.  If 
u{x,t)  =  (p{x  —  st)  is  to  satisfy  Oleinik’s  inequality  the  numerical  scheme  must  yield 
a  steady  progressing  profile  with  a  neirrow  transition  from  ur  to  ur  [21,  23].  Harten 
refers  to  this  property  as  resolution. 

If  the  solution  u{x,t)  =  (f[x  —  si)  is  inadmissable,  then  the  solution  is  a  fan  of 
waves  [23].  This  fan  of  waves  is  a  function  oi  xjt  and  consists  of  a  rarefaction,  or 
expansion,  wave  in  the  same  field  as  the  initial  discontinuity.  The  physical  solution 
requires  the  initial  discontinuity  break  up  instantaeously,  since  u{x,t)  —  <f){x  f t) .The 
term  entropy  enforcement  refers  to  the  requirement  that  the  numerical  scheme  break 
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up  the  initial  discontinuity  at  a  fast  rate,  thus  imitating  the  physical  behavior  [23] . 

The  systems  of  conservation  laws  under  consideration  contain  two  types  of 
characteristic  fields,  termed  nonlinear  and  linearly  degenerate  by  Harten  [23].  The 
nonhnear  fields  are  defined  such  that  a^R’‘  0,  while  the  linearly  degenerate  fields 

are  defined  by  =  0.  The  waves  of  a  nonlinear  field  are  shock  waves  or  expansion 
waves  while  the  waves  of  a  hneaxly  degenerate  field  are  solely  contact,  or  entropy, 
discontinuities. 

To  address  the  question  of  entropy  enforcement,  consider  the  scheme  given  by 
Eq  1.26,  which  has  the  effective  numerical  viscosity  coefficient 

/5K  A)  =  5  [«(>')->'’]  (!■«) 

The  least  dissipative  form  of  Q  is  arrived  at  by  choosing  it  to  be  consistent  with 
Eq  1.29  such  that 

(5(a;)  =  |x|  (1.50) 

With  Q  given  by  Eq  1.50,  the  scheme  of  Eq  1.26  can  be  rewritten  as  [23] 

=  u"  -  (uj+1/2)  -  (^Pj-1/2)  Aj_i/2u"  (1-51) 

where 

i/~  =  min(z^,  0)  =  ^  (1/ -  |i/|)  i/+ =  max(i/,  0)  =  |(i/ +  Iz^l)  (1-52) 

Harten  points  out  that  the  scheme  of  Eqs  1.51  and  1.52  is  a  generahzation  of  the 
Courant,  Isaacson,  and  Rees  scheme,  which  has  been  thoroughly  analyzed  in  the 
literature.  The  interested  reader  is  referred  to  reference  [23]  for  further  details. 

If  the  scheme  given  by  Eqs  1.51  and  1.52  is  apphed  to  the  Riemann  problem 
with  the  Rankine-Hugoniot  relation  satisfied  by  letting  the  speed  of  propagation  be 


15 


zero,  Eq  1.51  holds  the  initial  discontinuity  steady  regardless  of  entropy  consider¬ 
ations.  In  other  words,  the  initial  discontinuity  is  not  broken  up  and  there  is  no 
entropy  enforcement  in  this  case. 

The  problem  is  that  the  numerical  viscosity  vanishes  for  v  =  0.  Harten  elimi¬ 
nates  this  problem  by  modifying  Q{x)  =  |x|  near  x  =  0  to  be  positive.  The  modifi¬ 
cation  is  as  follows  [23]  for  0  <  e  <  | 


Q{x)  =  {x^/{4e))  +  e  |x|  <  2e 
=  |x|  |x|  >  2e 


(1.53) 


with  the  entropy  correction  parameter,  e,  typically  of  order  0.1. 

Harten  summarizes  the  results  of  numerical  experiments  carried  out  with  the 
scheme  of  Eqs  1.44  and  1.45  applied  to  the  Euler  equations  for  the  Riemann  problem. 
These  experiments  used  e  =  0.05,  0.1,  and  0.25  for  all  fields,  and  also  e  =  0  for  the 
linearly  degenerate  field.  Basically,  highly  resolved  shocks  were  obtained  for  all  values 
of  e  under  consideration.  The  contact  surface  was  better  resolved  than  with  the  first- 
order  accurate  scheme  of  Eqs  1.51  and  1.52,  but  stiU  remained  rather  smeared. 

To  prevent  excessive  smearing  in  the  linearly  degenerate  field  containing  the 
contact  surface,  Harten  replaces  Eq  1.46  in  the  hnearly  degenerate  field  with  [23] 

9j  =  9j  +  Sj9j  (1-54) 

where  gj  is  the  right  hand  side  of  gj  given  by  Eq  1.46  and,  gj  is 


gj  =  S max  [o,min  o’i+i/2  Q:j+i/2|)] 

(1.55) 

S  =  sgn  (ay+i/z) 

(1.56) 

o-j+i/2  —  o-i+i/2  (^'2+1/2)  ~  ^  ('^i+1/2)] 

(1.67) 
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ai+1/2 


-“i-i/2  /  I  «i+i/2 


+ 


«i-l/2|) 


(1.58) 
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11.  Analysis 


2.1  Euler  Equations 

The  Euler  equations  are  statements  of  the  conservation  laws  for  mass,  momen¬ 
tum,  and  energy  assuming  an  inviscid,  nonconducting  gas.  When  the  Euler  equations 
are  arranged  such  that  p,  pu,  pv,  and  e  cire  the  dependent  variables  the  conservative 
or  divergence  form  is  obtained.  Leix  showed  that  the  conservative  form  of  the  Eu¬ 
ler  equations  satisfies  the  weak  solution  of  the  Rankine-Hugoniot  relations  and  thus 
correctly  predicts  the  jump  conditions  across  the  shock  discontinuity  [1,  34].  In  fact, 
use  of  the  conservative  form  is  necessary  for  the  discontinuity  to  represent  a  physical 
wave  when  shock  capturing  schemes  are  applied  [1].  The  conservative  form  is  often 
referred  to  as  the  divergence  form  because  the  equations  identify  the  divergence  of 
physical  quantities.  The  governing  equations  may  be  written  in  the  following  vector 
form; 

where  U  contains  the  dependent  VEiriables  which  are  the  density,  p\  s-momentum, 
pu]  y-momentum,  pv\  and  total  energy  per  unit  volume,  e.  F  contains  the  flux  terms 
differentiated  with  respect  to  x,  and  G  contains  the  flux  terms  differentiated  with 
respect  to  y.  The  elements  of  t/,  F,  and  G  are: 


p 

m 

n 

m 

F  = 

rrF  Ip  +  p 

G  = 

nu 

n 

mv 

Ip  -\-p 

e 

(e  -1-  p)m/p 

.  (e  +  p)’^/p  . 

where  m  =  pu  and  n  =  pv.  The  pressure,  p,  is  given  as 


p  =  (7  -  1) 


(m2  +  n2) 

2p 


(2.2) 


(2.3) 
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for  a  thermally  and  calorically  perfect  gas. 

A  general  spatial  transformation  of  the  form  ^  =  C{^,y)  and  rj  —  r]{x,y)  is  used 
to  transform  Eq  2.1  from  the  physical  domain  {x,y)  to  the  computational  domain 


(^,77).  The  strong  conservation  law  form  of  the  Euler  equations  is 

now  given  by  [42] 

dir  aHu)  da{u) 

dt  di  dr) 

(2.4) 

U^UIJ 

(2.5) 

F  =  ((.F  +  4G)  jJ 

(2.6) 

a  =  {r,.F  +  ,,yG)/J 

(2.7) 

J  —  ^yVx 

(2.8) 

where  J  is  the  Jacobian  of  the  transformation. 

Since  the  TVD  method  used  herein  utilizes  the  local-characteristic  approach, 
which  is  a  generahzation  of  Roe’s  approximate  Riemann  solver[35],  the  Jacobians  A 
and  B  of  F  and  G  are  required  and  can  be  written  as 


(2.9) 

B  =  (j),4  +  nyB) 

(2.10) 

where 


A  =  Fu  = 


0  1 

1(7  -  1)  {u^  +  v^)  -u^  (3  -  7)u 

—uv  V 


0  0 

(1  —  j)v  7  —  1 
u  0 

(1  —  'y)uv  7^ 
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0 


0 


1 


0 


B  —  Gu  = 


—uv 

|(7  -  1)  (u^  +  v^)  - 


V  u  0 

(1  — 7)u  {3  — 'y)v  7  — 1 

(1  —  'y)uv  H  —  {'y  —  l)v'^  jv 


with  the  total  enthalpy,  H,  given  by 


H  = 


IP 

(7  -  l)/> 


The  eigenvalues  of  A,  denoted  {a}^,  ,  are 


^  ixU  +  iy'v  -  ^ 

$xu  +  CyV 

ixU  +  ^yV  +  k^c 
\  (xU  +  ^yV  / 


(2.11) 


(2.12) 


where 

k  =  + 

The  right  eigenvectors  of  A,  i?|,  Rfj,  are 


1 

1 

u  —  kic 

R\  = 

u 

V  —  k2C 

V 

H  —  k\uc  —  k2vc 

\{u^  +  V^)  _ 

1 

0 

u  +  kiC 

II 

-k2 

V  +  k2C 

h 

H  +  kiuc  +  k2vc 

k\v  —  k2U 

(2.13) 


(2.14) 
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where 


and 


(2.15) 


(2.16) 


The  eigenvalues  and  eigenvectors  of  B  are  obtained  by  replacing  ^  in  Eqs  2.12 
through  2.16  with  tj: 


^  TJxU  +  TjyV  -  kr,c  '' 
Tf^U  +  T]yV 
T]^U  +  r]yV  +  kr,c 


\  ilxU  +  rjy'^  J 


11 

l  +  vl 

1.0 

1.0 

u  —  kiC 

u 

V  —  k2C 

V 

H  —  kiuc  —  k2vc 

_  l{u^  +  v^) 

1.0 

0.0 

u  +  kic 

K  = 

-k2 

V  +  k2C 

h 

H  +  kiuc  +  k2vc 

kiv  —  k2U 

h 

k2 


(2.17) 


(2.18) 


(2.19) 


(2.20) 

(2.21) 
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2.2  Numerical  Procedure 


2.2.1  2-D  Harten-Yee  Finite- Volume  Algorithm.  An  upwind  TVD  scheme 

in  finite-volume  form  [42]  is  used  in  the  present  study.  The  grid  spacing  is  denoted 
by  and  At],  such  that  ^  —  jA^  and  rj  =  kArj.  Utihzing  the  Strang-type  fractional 
step  method  allows  the  scheme  to  be  implemented  in  a  local-characteristic  approach 
and  ensures  second-order  accuracy: 


where 


(2.22) 

=  (>h  =  -  Niji) 

(2.23) 

(2.24) 

with  h  =  At.  Apphcation  of  the  entire  sequence  of  operators  (one  iteration)  advances 
the  solution  two  time  levels.  The  functions  ^  and  ^  numerical  fluxes 

in  the  (  and  rj  directions  evaluated  at  cell  interfaces.  For  instance,  Fj+i,k  Yee’s 
finite-volume  formulation  is  expressed  as 


^f)  +  Gj+i,k) 


(2.25) 


where  the  subscript  ^  -1-  |  is  a  simpHfied  notation  for  j  |,fc.  The  numerical  flux 
function  in  the  t]  direction  is  defined  similarly: 


(2.26) 


The  eigenvalues  and  eigenvectors  are  evaluated  at  cell  interfaces  using  sym¬ 
metric  averages  of  Uj^k  and  Uj+i^k,  Uj^k  and  Uj^k+i,  respectively.  Roe’s  averaging 
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technique  for  a  perfect  gas  is  used  herein  and  takes  the  form  [35] 


Duj^l^k  '^j,k 


Ht  •  I  1  ^  __ 

^+2’*=  D  +  1 


Dvj+i^k  +  Vj^k 
-  D  +  1 


DHj^i^k  +  Hj^k 


D  +  1 

\,k  =  (7  -  1)  -  2  (“?+i,fc  + 


(2.27) 

(2.28) 

(2.29) 

(2.30) 

D  =  \l Pi -^i,kl  Pi, k  (2.31) 

Roe’s  averaging  technique  is  used  because  it  has  the  computational  advantage  of 
perfectly  resolving  stationary  [43],  but  not  necessarily  moving,  discontinuities. 

The  quantities  and  l/^y+i  are  defined  as  follows  for  the  finite-volume 

formulation 


where 


'i+k 


1 

2 


'1^  ^ 


/  i+l,k 


1  _  1  /  1  ^  1 
^i+k  2  Ji,k, 


(2.32) 

(2.33) 


The  constants  and  {k2)j^  necessary  in  determining  Eq  2.15,  are 

defined  as 


i+l 


\/(j)i+i  +  (j),- 


(2.34) 


i+ 


and 


(fe)i 


i+ 


'(^)L + (^)m 


(2.35) 
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The  vector  function  is  composed  of  elements  denoted  as 

second-order  upwind  TVD  scheme.  The  elements  are  given  by 


^  (^5+^  +  4)  -  <5  (^'+1  +  7j+t)  « i+L 

where,  with  A  =  At/A^ 

•'iH  = 

a-,  I  is  the  difference  of  the  characteristic  Vciriables  in  the  £  direction, 

3-r  2 


aj+t  =  R-li  iUj+i,k  -  Uj^k) 


(2.36) 


(2.37) 


(2.38) 


or 


3 

i+J 

^  4 

.  i+|  . 

L 

where 


aa 


7 


r-2 

i+| 


(aa  —  bb)/2 
Aj^p  -  aa 
(aa  bb)  1 2 
cc 


A  .+  16  + - ^A.+ip  -  u-^A.^m  -  ^,-+iA  .+  Ln 


(2.39) 


(2.40) 


bb  =  — ^  ^kiAj^im  -  (kiUj^i  +  ^z'Wj+i)  Aj+ip  -f-  k2Aj^in^  (2-41) 

cc  =  kiAj_^in  +  {k2u-j^i_  -  fcivj+i)  Aj.+  lP  -  fc2Aj-+Lm  (2.42) 

Aj.+  p  =  Zj+i,k  -  Zj,k  (2.43) 

The  difference  of  the  local  characteristic  variables  in  the  rj  direction  is  obtained  in 
similar  fashion: 

“i+l  =  JliU  (UjMi  -  Ui.k)  (2.44) 


24 


((id  —  ee)/2 

“fc+J 

Ajt+ip  —  aa 

(dd  +  ee)/2 

where 


dd  = 


7-1 


Afc+tp  -  -Ufc+iAfc+im  -  Ufc+iAjt+Ln 


ee 


with 


=  [^iA;i+im  -  (fciMfc+i  +  ^2’^fc+L)  Afc+ip  +  fe2A;.+in] 

//  =  ^1  Afc+in  +  (A:2«fc+i  -  hvk+i)  Afc+i/J  -  A:2A;t+Lm 


(^l)fc+l  - 


(^2)fc+L  = 


^2*V  4.  ( U3L 

Jk+ 


k+k 


i^)k,, 


A(^_l_lZ  —  ^j,fc+l  Zj^k 


The  function  7'  ^  is  given  by 

Ji"  3 


=  < 


KO^ 


“i+l 


where 


o-(®)  =  ^  [<5(®)  - 


(2.45) 


(2.46) 

(2.47) 

(2.48) 

(2.49) 


(2.50) 


(2.51) 


(2.52) 


(2.53) 
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and 


Q(x}  = 


\x\ 


(2.54) 


(|i|  >  2e) 

(*7(4e))  +  e  (li!|<2e) 

The  entropy  correction  parameter,  e,  is  generally  fixed  during  computations,  but  can 
vary  between  0  and  0.5. 


The  function  g^-  in  Eq  2.36,  initicilly  referred  to  in  Section  1.2.2,  is  termed  the 
‘hmiter’  function  and  can  be  expressed  in  a  variety  of  ways  [42].  The  present  study 
bases  the  choice  of  the  Hmiter  on  the  type  of  characteristic  field  under  consideration. 
For  the  nonhnear  fields,  a}E}  ^  0,  Eq  4.34d  of  Yee  [42]  is  used: 


+ 


(2.55) 


For  the  Hnearly  degenerate  fields,  =  0,  Eq  4.34g  of  Yee  [42]  is  appHed: 


=  S  ■  max  0,  min  ^2  ,  S  ■  >  min  ,  2S  ■  Q:j_L) 


(2.56) 


where 

S  =  sgn  (a'^i) 


(2.67) 


The  nonhnear  fields  correspond  to  /  =  1  and  1  =  3  while  the  Hnearly  degenerate 
fields  correspond  to  /  =  2  and  1  =  4.  It  should  again  be  noted  that  the  waves  of  a 
nonlinear  field  are  either  shocks  or  raxefraction  waves  while  the  waves  of  a  Hnearly 
degenerate  field  are  uniquely  contact  discontinuities  [23].  Since  this  is  a  five-point 
scheme,  the  values  of  are  needed  at  cell  centers  just  outside  the  computational 
domain.  Zeroth-order  extrapolation  is  used  to  obtain  the  necessary  values,  foUowing 
the  example  of  Harten  [23]. 


2.2.2  2-D  Harten-Yee  Chain-Rule  Algorithm.  In  addition  to  the  finite- 

volume  formulation  of  Yee,  a  finite-difference  form  based  on  the  chain-rule  con- 
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servation  form  of  the  governing  equations  was  utilized  [38]: 


m  dF{u) 

dt  di 


+  ^y — - ’7= 


di 


dr] 


dG{U) 

drj 


=  0 


(2.58) 


Previous  researchers  report  that  the  governing  equations  in  this  form  are  more  com¬ 
putationally  efficient  than  the  strong  conservation  form  used  in  the  finite-volume 
approach  [38].  This  was  found  not  to  be  the  case  for  the  current  TVD  algorithms, 
with  both  formulations  performing  approximately  the  same  in  terms  of  computa- 
tioncd  efficiency. 


The  local  characteristic  approach  given  by  Eq  2.22  is  now  applied  to  U  instead 
of 


(2.59) 

where 

=  U’,  = 

(2.60) 

(2.61) 

The  numerical  fluxes,  and  ,  for  the  chain-rule  conservation  form  are 

2 

{Cx)j,k  {Fj,k  +  T'i+i.fc)  +  (4)i,fc  (G'i.fc  +  Gj+i,k)  + 

(2.62) 

and 

~  2 

An 

{Vx)ik  (Rj.k  +  Fj^k+i)  +  (J7y)j,fc  {Gj.k  +  Gj.fc+i)  + 

(2.63) 

The  quantities  and  (^2)^+1  are  defined  as  follows  for  the  chain  rule 

formulation: 

=  J  [(e.)i,.  +  (2.64) 
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=  - 

(2.65) 

=  - 

(^y)i+l 

(2.66) 

2.3  Navier-Stokes  Equations 

The  conservative  form  of  the  Navier-Stokes  equations  is  written  as 


^  dF{U)  dG{U)  dF,{U,U,,Uy)  dG4U,U.,Uy) 

dt  dx  dy  dx  dy 


where  U,  F,  and  G  are  the  same  as  for  the  Euler  equations,  Eq  2.1.  and  G„  are 
the  viscous  flux  terms,  given  as 


0 

0 

Txx 

II 

Txy 

'^xy 

"^yy 

UTg,*  +  -  qx  _ 

_  UTxy  -f  VTyy  -  qy  _ 

(2.68) 


Txx,  7'xii)  and  Tyy  are  the  viscous  stresses: 


Txx  —  “i” 

Txy  =  n{uy  +  v^)  (2.69) 

Tyy  —  (2^  “t"  \^Vy  “j"  AUg; 


where  y,  and  A  are  the  first  and  second  coefficients  of  viscosity  respectively.  The  first 
coefficient  of  viscosity  is  determined  using  Sutherland’s  formula  [1]; 


/x  =  Gi 


J<3/2 

T  +  C2 


(2.70) 
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where  Ci  -  1.458  x  10“®  kgj  {m-  s-  \/^)  and  C2  =  110.4  K.  The  second  coefficient 
of  viscosity  is  given  by 

B-2  +  -  (2.71) 

where  B  =  Aj 3  yields  Stoke’s  hypothesis,  A  =  -2/3/z.  Solutions  were  also  arrived 
at  using  B  =  2,  based  on  Sherman’s  work  as  reported  by  White  [45].  No  difference 
was  observed  in  the  numerical  solutions  using  5  =  4/3  or  5  =  2. 


The  quantities  and  are  components  of  the  heat  flux  vector,  q  =  —kVT. 
The  coefficient  of  thermal  conductivity,  k,  is  determined  from  the  Prandtl  number, 
Pr: 


Pr  = 


k 


(2.72) 


with  Pr  =  0.72  for  air. 


The  equations  may  be  written  in  Hnearized  form  as 


Ut  +  AUx  +  BUy  =  AlUx  +  BlUy  +  A2UXX  +  B2Uyy  +  (^3  +  B3)  U xy  (2.73) 

where  the  viscous  Jacobian  matrices  are 

A^^dF.IdU  A2^dFjdUx  A3==dFjdUy 

B,  -  dGJdU  B2  =  dGJdUy  B3  =  dGJdUx 

with  the  individual  terms  given  in  Appendix  A. 

A  general  spatial  transformation  of  the  form  ^  and  77  =  T][x,y)  is  used 

to  transform  Eq  2.73  from  the  physical  domain  (s,j/)  to  the  computational  domain 

{LnY- 

Ut  +  AU^  +  BUr,  =  AiU(  +  BiUr,  +  A2U^^  +  B2Ur,rj  +  (^3  +  BaJ  U^r)  (2.75) 
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where 


and 


A  =  + 

A-l  =  +  CyBi 

A-2  =  CIA2  +  ^lB2  +  (^3  +  B3) 

A.3  =  ixIlyAz  +  CyVxBs  +  ^xVxA2  +  ^y^y-^2 

B  =  tj^A  +  rjyB 

Bl  =  -qxAl  +  TJyBl 

B2  =  vl^2  +  VIB2  +  VxVy  +  B3) 

B3  =  ^xVyB3  +  ^y?7x^3  +  ^x77x^2  +  ^y?7y^2 


(2.76) 


(2.77) 


2.4  Numerical  Procedure 

A  first-order  time,  second-order  space,  upwind  TVD  scheme  is  now  presented 
for  the  Navier-Stokes  equations.  Based  upon  the  excellent  results  achieved  in  the 
inviscid  case,  a  chain-rule  formulation  is  utihzed.  The  scheme  for  the  Euler  equations 
,  described  in  Section  2.2.2,  is  second-order  accurate  in  space  and  time.  Taylor  series 
expansion  shows  the  scheme  is  a  representation  of 


Ut  +  ^xF(  +  T]:,Fr,  +  ^yG^  +  r]yG,j  =  ^^-Utt  + A'^U(^+ (^AB  +  BAjU^r, 

+B"£/„]+0[A^^A^^A,“) 


(2.78) 


and  is  second-order  accurate  for  the  Euler  equations,  since 


Utt  =  A^Ui^  +  {Ab  +  bA)  Uirj  +  B^Ur^r,  (2.79) 


Viscous  terms  are  added  to  the  Euler  scheme,  Eqs  2.60  and  2.61,  using  second- 
order  accurate,  central-difference  approximations: 


At 


{F;^,,-F;_.J+At^-, 


(2.80) 


30 


At 
“  Arj 


ctuu  =  Uh  -  ~  (o-^L  -  +  Ai®; 


-  '^i.k  ~  ~  ^  (2-81) 

where  F  and  G  are  given  by  Eqs  2.62  and  2.63.  The  viscous  terms,  and  are 


5^=2!^  (ft,,..  -  ft,,.)  +  (f,),.  (G.,„,.  -  G,,,)]  (2.82) 


and 


(2.83) 


In  order  to  maintain  second-order  accuracy,  the  derivatives  appearing  in  the 
viscous  terms  of  and  G„  must  be  differenced  appropriately.  The  C  derivative 
terms  in  F„  and  C?„  are  backward  differenced  when  computing  Similarly,  the  77 
derivative  terms  in  F„  and  G„  are  backward  differenced  when  computing  AU 
other  derivative  terms  are  central  differenced. 


The  scheme  given  by  Eqs  2.80  through  2.83  is  a  representation  of 


Ut  +  ^xF(  -f  T]xFrj  -f  +  VyGr)  —  -f  -f  CyG^^  -f  T^yG^r, 

-t-f  [-Utt  +  A^Ua  +  [aB  +  BA) 
+B^Ur„]+0[At\Ae,A7i^] 

(2.84) 


Examination  of  Eq  2.75  reveals 


Uu  =  {B^  +  Bl-BB,-BrB)Ur,r, 

^AiBi  -f  BiAi  —  A\B  —  bAi  —  ABi  —  B\A  -f  Ab  -i-  bAl^  u^ri 

+  {B\B2  -I-  B2B1  —  BB2  —  B2B^  ^■nriT)  +  ^2^rinn'n 
-j-  (^BsBi  -f-  B1B3  —  B3B  —  BB3  -j-  B1A3  A3B1 
— BAs  —  A3B  B2A1  -|-  A1B2  —  B2A  —  AB^  ^inn 
-f  [B2BZ  -\r  B3B2  -{■  B2A3  -f  A3B-^ 

+  -|-  Ag  A2B2  -(-  B2A2  -|-  A3B3  4- 
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(2.85) 


+  4"  A\A.3  —  -AzA  —  AAz  “h  AiBz  4"  B^Ai 

— ABz  —  Bz-A  4"  AzBi  4'  B-^Az  —  AzB  —  BA^ 

+  ^.4.2.43  4"  AzAz  4-  ^42.63  4“  BzAzJ 

+  (i=*  +  i?-iii-ixi)f/« 

4-  {^AiAz  4-  AzAi  —  AAz  —  -42.4^  4- 

Since  the  term  of  0[At]  in  the  truncation  error  of  Eq  2.84  does  not  vanish  upon 
substitution  of  Eq  2.85,  the  scheme  obtained  by  adding  central-difference  represen¬ 
tations  of  the  viscous  terms  to  the  Euler  scheme,  Eqs  2.80  through  2.83,  is  first-order 
accurate  in  time  and  second-order  accurate  in  space.  This  new  scheme  represents 

Ut  4-  4-  tIxFt}  4-  $yG(  4-  %G,j  =  ^xFv^  4-  V^Pv,,  4-  CyGv^  4-  ^lyGv^ 

+0[At,Ae,^ri^] 

This  first-order  time,  second-order  space  scheme  is  hereafter  referred  to  as  ATNSCl 
which  is  shorthand  for  Ist-Order  AFIT  TVD  Navier-Stokes  Code. 


2.5  Baldwin- Lomax  Turbulence  Model 


Since  this  effort  represents  a  first  step  in  determining  the  performance  of  TVD 
schemes  in  conjuction  with  a  turbulence  model,  the  Baldwin-Lomax  algebraic  turbu¬ 
lence  model  [3]  was  chosen  for  simplicity.  The  effects  of  turbulence  are  simulate  by 
adding  an  eddy  viscosity  fit  to  the  molecular  viscosity  fi,  with  the  total  viscosity  thus 
given  by  /i  -t-  /if.  The  heat  flux  terms  Eire  cailulated  by  replacing  k/Cp  =  fi/Pr  with 
k/Cp  =  fi/Pr  4-  fitlPrt-  It  is  a  two  layer  model  with  the  turbulent  eddy  viscosity 
given  by 


(^t) 


inner 


outer 


(y  —  y  crossover) 
{y crossover  ^  V) 


(2.87) 


where  y  is  the  distance  away  from  the  surface,  in  the  normal  direction.  The 
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value  of  y  at  which  the  inner  and  outer  eddy  viscosites  are  first  equal  is  designated 

y crossover 


The  Baldwin-Lomax  model  utilizes  the  following  formulation  for  the  inner  re¬ 
gion  eddy  viscosity 


(/^Oinner  = 


(2.88) 


with  the  length  scale  given  by 


1  =  ky^l-  exp 


(2.89) 


the  vorticity  calculated  as 


m  = 


du  dv 


\dy  dx\ 

and  the  boundary  layer  coordinate  y'^  defined  as 


(2.90) 


^  ^/Pw'^wV 
= - 

Pw 


The  outer  region  eddy  viscosity  is  calculated  according  to 


(2.91) 


{Pt)  outer  “  ^(^cpPFyjakeFKlehiy) 


(2.92) 


with 


F,u^ke  =  min 


ymaxF max 

C why max"^ dif  I Frr 


Fmax  is  the  maximum  of  the  function 


(2.93) 
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F{y)  =  ykl  [l  -  exp  (-y+/A+)] 


(2.94) 


and  ymax  is  the  value  of  y  at  the  location  where  F[y)  is  a  maximum.  The 
exponential  term  of  Eq.  2.94  is  taken  to  be  zero  in  a  wake.  FKUb  is  the  Klebanoff 
intermittency  factor 


FKUbiy)  =  [l  +  5.5  {CKuiy/ymaxf] (2.95) 

The  difference  between  the  mciximum  and  minimum  velocity  in  a  given  profile, 
at  a  fixed  x  location,  is  termed  Udif  and  is  given  by 


Uiif  =  [s/u^  +  vA  —  [s/u^  +  vA  ,  (2.96) 

The  minimum  value  of  velocity  is  zero  when  the  profile  begins  or  ends  at  a 
wall, and  is  therefore  nonzero  only  in  wakes. 

Two  different  transition  criteria  are  used  in  the  current  study.  The  Baldwin- 
Lomax  criteria  is  to  set  the  eddy  viscosity  equal  to  zero  when  the  following  condition 
is  met; 

(/Xt)  <  Cjnuimy'oo  (2.97) 

The  second  criteria  is  termed  the  zero  shear  stress  test  and  is  used  based  on 
Boyle’s  [4]  reccomendation.  Simply  stated,  transition  is  assumed  to  occur  when  the 
shear  stress  reaches  zero  at  a  surface. 

The  constants  used  in  the  Baldwin-Lomax  model  aire  [3] 
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A+ 

=  26 

C'.p 

=  1.6 

CKlth 

=  0.3 

C^k 

=  0.25 

k 

=  0.4 

K 

=  0.0168 

Pvt 

=  0.9 

P'  muim 

=  14 

(2.98) 
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III.  Results  and  Conclusions 


3.1  Boundary  Conditions  for  the  Inviscid  Studies 

Appropriate  boundary  conditions,  in  conjunction  with  initial  conditions  and 
flow  parameters  such  as  Mach  number,  are  necessary  to  arrive  at  the  particular 
solution  of  interest.  The  boundary  conditions  are  now  described  in  detail. 

3.1.1  Inlet  and  Exit  Boundary  Conditions.  If  the  inlet  velocity  is  supersonic, 
aU  characteristics  originate  upstream  of  the  computational  boundary  so  the  four  nec¬ 
essary  flow  quantities  may  be  specified.  Likewise,  if  the  outflow  velocity  is  supersonic 
aU  characteristics  originate  inside  the  computational  domain  and  the  four  necessary 
exit  quantities  must  be  extrapolated  from  the  interior.  Second-order  accurate  ex¬ 
trapolation  is  utiUzed  in  the  schemes  under  consideration. 

Subsonic  inflow  and/or  outflow  presents  a  more  compHcated  situation.  In  ap¬ 
plying  the  boundary  conditions  at  the  inlet  and  exit  of  the  domain,  it  is  assumed 
that  these  boundaries  are  sufficiently  distant  from  the  cascade  so  that  planar  wave 
disturbances  propagate  coUinearly  with  the  stream  function.  The  disturbances  are 
required  to  leave  the  computationcd  domain  without  reflection,  except  for  the  re¬ 
flection  of  pressure  disturbances  at  the  exit.  For  subsonic  inlet  velocities,  the  inlet 
boundary  conditions  are  arrived  at  by  first  assuming  that  the  inlet  is  part  of  an 
imaginary  duct  extending  infinitely  far  upstream  of  the  cascade.  AU  waves  radiating 
from  the  computational  domain  should  pass  the  inlet,  without  reflection,  and  con¬ 
tinue  travelling  upstream  for  aU  time.  Specification  of  a  constant  thermodynamic 
state  at  upstream  infinity  requires  the  expansion  disturbance  traveUing  upstream  to 
behave  as  a  simple  wave.  This  behavior  aUows  the  application  of  one-dimensional 
characteristic  theory  at  the  inlet  [17]. 

For  subsonic  inflow,  only  one  characteristic  runs  from  the  interior  of  the  domain 
towards  the  computational  boundary.  Therefore,  three  quantities  must  be  specified 
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while  one  may  be  extrapolated  from  the  domain  interior.  Fcir  upstream,  the  total 
pressure,  and  toted  temperature,  Ti^,  are  specified,  while  only  the  inlet  flow 
angle,  ^2,  is  specified  at  the  computational  boundary.  The  speed  of  sound  at  the 
inlet,  C2,  is  extrapolated  from  the  domain  interior.  The  Riemann  invariant  along  the 
chciracteristic  spanning  the  expansion  wave  from  leading  to  trailing  edge  is  given  by 

Voa  H - -Coo  =  14  H - rC2  (3.1) 

7-1  7-1 

where  V  is  the  magnitude  of  the  velocity  vector.  As  the  velocity  vanishes  far  up¬ 
stream,  the  inlet  velocity  is  obtained  from 

^  {c„  -  C,)  (3.2) 

7-1 

which,  along  with  the  inlet  flow  angle,  determines  u  and  v.  The  inlet  pressure  is 
determined  from  the  isentropic  relation 

—  (3.3) 

The  speed  of  sound  and  pressure  fix  the  state  point,  uniquely  determining  the  density 
and  internal  energy. 

For  subsonic  axial  Mach  numbers,  simple-wave  theory  is  also  apphed  at  the 
exit.  The  exit  is  treated  as  an  open-end  duct  that  exhausts  into  a  plenum,  requiring 
the  exit  pressure  to  match  the  plenum  pressure.  Thus,  aU  pressure  disturbances 
are  reflected  back  into  the  computationed  domain  from  the  exit.  Two  characteristics 
extend  from  the  interior  of  the  computational  domain  to  the  exit,  while  one  originates 
outside  the  domain.  Thus  only  one  quantity,  in  this  case  pressure,  can  be  specified  at 
the  exit.  AU  other  quantities  must  be  extrapolated  from  the  interior  of  the  domain. 
The  quantities  chosen  for  extrapolation  are  entropy,  tangential  velocity,  and  the 
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Riemann  invariant,  The  density  is  obtained  from  the  isentropic  relation 


Pz  = 


(3.4) 


where  Sint  is  the  entropy  extrapolated  from  the  interior.  The  pressure  and  density  fix 
the  state  point,  uniquely  determining  the  speed  of  sound  and  internal  energy.  With 
the  tangential  velocity  extrapolated  from  the  interior,  the  axial  velocity  is  obtained 
by  applying  the  Riemann  invariant  in  the  axial  direction: 

2 

U3  =  Rinv - -C3  (3-5) 

7-1 


where 

2 

Rinv  —  T  Qnt 

7  -  1 

and  Uint  and  c,-„t  are  the  axial  velocity  and  speed  of  sound  at  the  point  inside  the 
domain  where  the  Riemann  invariant  is  evaluated. 


3.1.2  Periodicity  and  Blade  Boundary  Conditions.  Only  one  blade  passage  of 
an  infinite  cascade  is  analyzed.  Therefore,  periodicity  conditions  are  apphed  at  cell 
centers,  or  ghost  points,  located  outside  the  computational  domain.  These  points 
are  located  cdong  the  outer  boundary  and  also  along  the  wake  cut  when  a  C-type 
grid  is  utiHzed.  For  an  H-type  grid,  ghost  points  are  located  along  the  upper  and 
lower  boundaries  upstream  and  downstream  of  the  blade.  At  the  blade  surface,  the 
only  condition  that  can  be  specified  is  the  requirement  for  surface  tangency.  Since 
the  blade  surface  is  mapped  to  a  constant  77  coordinate,  the  normal  component  of 
velocity  is  given  by 


rjxU  +  'TJyV 


while  the  tangential  component  is 


y  ^ 

\/vl  +  nl 

The  requirement  for  surface  tangency  is  met  by  setting 


(3.8) 


(3.9) 


and 

Ki,.  =  -K.,-,  (3.10) 

where  j  is  the  ^  index,  0  represents  a  ghost  point  just  inside  the  body,  and  1  is  the 
index  of  the  first  cell  center  above  the  body.  Cell  centers  and  ghost  points  are  used  to 
place  the  blade  surface  along  the  interface  of  the  grid  cell  and  ghost  cell.  This  mesh 
system  helps  ensure  both  consistent  and  conservative  boundary  conditions  [34].  The 
inverse  relation  between  the  Cartesian  velocities  and  Eqs  3.8  and  3.7  then  gives 


Uj,0 

1 

Vy 

Px 

1 

0 

i^j,0 

Py  . 

1 

0 

(3.11) 


The  pressure  at  the  ghost  points  is  obtained  by  applying  the  normal-momentum 
equation  at  the  first  fine  of  cell  centers  above  the  body  [32]: 

-p  +  lyu)  {t]^u^  +  r]yV^)  =  (77^:6  +  iyily)  Pi  +  ifi  +  ^y)  Pn 

=  Pn^Jpl  +  (3.12) 

Central  differences  are  used  for  both  the  ^  and  77  derivatives. 
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One  additional  property  is  needed  to  fix  the  state  of  the  ghost  points.  In  the 
present  study,  an  adiabatic  wall  condition  is  chosen  to  provide  this  final  property: 

=  0  (3.13) 

Although  it  is  inconsistent  with  the  Euler  equations  to  specify  either  the  temperature 
or  its  gradient  at  the  blade  surface  [34],  an  adiabatic  waU  condition  has  been  used 
by  others  [32]  and  yields  results  that  agree  weU  with  theory  and  experiment. 

3.2  Boundary  Conditions  for  Viscous  Flow 

Boundary  conditions  for  viscous  flows  are,  in  general,  more  straight  forward 
than  their  inviscid  counterparts  described  in  Section  3.1.  At  the  wall,  the  invis- 
cid  surface  tangency  condition,  Eqs  3.7  and  3.8,  is  replace  by  the  viscous  no-sHp 
requirement: 

(3.14) 

u  =  0 

Simphfied  wall  temperature  contitions  representing  either  an  adiabatic  wall 

=  0  (3.15) 

or  constant  temperature  waU 

T^  =  C  (3.16) 

are  used  in  the  current  study,  depending  on  the  flow  of  interest.  With  the  waU 
mapped  to  a  constant  rj  coordinate,  the  pressure  at  the  waU  is  obtained  by  solving 
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the  normal-momentu.m  equation: 


Pn^jril  +  Vl  =  i^xrjx  +  ^yr]y) P(  +  {vl  +  Vl)  Pr) 

=  +  Py'^n)  {^»  +  ’nxK) 

+^y  [4  (2^  +  +  Vy  +  •^)  J  }  ^3 

“t"  Vy'^v)  {^®  i^yPC  "t"  VyPn)  “I"  Py  V^^P’i)^ 

+PVx  +  "^^yVyUiv  +  Pl^rtr,) 

+  {2p  +  X)T]y  [dvii  +  KVyf^ir,  +  vlvrin) 

Flow  at  the  inlet  and  exit  of  the  computational  domain  is  assumed  to  be 
inviscid.  Inflow  and  outflow  relations  from  Section  3.1  are  thus  used  to  determine 
flow  quantities  at  these  boundaries.  As  stated  in  Section  3.1,  for  supersonic  outflow 
aU  quantities  must  be  extrapolated  from  the  interior  of  the  domain.  In  practice, 
this  extrapolation  is  also  performed  in  the  subsonic  boundary-layer  embedded  in  the 
supersonic  outflow.  For  the  cases  to  be  considered  herein,  no  adverse  effects  of  this 
extrapolation  are  noted. 

3.3  Shock-Boundary  Layer  Interaction 

An  indepth  experiment  in  laminar  shock-boundary-layer  interaction  was  car¬ 
ried  out  by  Hakkinen  et  al.  [19]  in  1959  at  the  Massachusetts  Institute  of  Technology 
under  the  sponsorship  of  the  National  Advisory  Committee  for  Aeronautics.  De¬ 
tailed  measurements  were  made  of  pressure  distribution,  skin  friction  coefficient, 
and  velocity  profiles  for  a  number  of  combinations  of  overall  pressure  ratio,  Pfipoo, 
and  shock  Reynolds  number,  Re^,,  at  a  freestream  Mach  number  of  2.0  for  a  shock 
wave  impinging  upon  a  flat  plate  boundary-layer. 

The  experimented  pressure  and  skin  friction  profiles  for  this  case  are  shown  in 
Fiqure  3.1,  and  a  sketch  of  the  wave  structure  is  shown  in  Figure  3.2.  The  friction 
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coefficient,  (7/,  is  defined  as 


(3.18) 


where  is  the  normcd  component  of  sheeir  stress  at  the  wall 

(3.19) 

and  Qoo  is  the  dynamic  pressure 

fe  =  (3.20) 

With  the  tangential  velocity  given  by  Eq  3.8,  and  the  wall  mapped  to  arj  =  constant 
coordinate,  can  be  written  as 

li{VyUr,-r]xVr,)  (3.21) 

No  negative  values  of  skin  friction  are  shown  because  the  total-head  tube  was  not 
able  to  reUably  indicate  negative  shear  values. 

A  typical  grid  used  for  the  numerical  investigations  is  shown  in  Figure  3.3. 
Spacing  is  held  constant  in  the  axial  direction  at  Axjxfhock  =  0.013  and  the  minimum 
spacing  in  the  normal  direction  is  dicated  by  specifing  =  2.2  per  Boyle  [4].  The 
computational  domain  is  initialized  at  the  uniform  freestream  conditions  to  the  left 
of  the  point  along  the  upper  boundary  at  which  the  shock  is  generated.  Post-shock 
conditions  me  apphed  downstream  of  this  point.  An  adiabatic  wall  condition  is  used 
to  obtain  the  wall  temperature  along  the  plate  and  the  nomal  momentum  equation  is 
solved  to  obtain  the  wall  pressure  in  combination  with  the  no-shp  velocity  constraint 
at  the  wall. 

Figure  3.4  depicts  the  solution  obtained  with  e  =  0  in  the  nonlinear  fields  and 
e  =  0.025  in  the  Hnearly  degenerate  fields,  based  on  the  author’s  previous  work  /cit- 
edr:91.  The  pressure  profile  of  Figure  3.4  clearly  shows  the  pressure  rise  to  separation. 
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Figure  3.2.  Flowfield  Structure 
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Figure  3.3.  Grid  Used  in  Shock-Boundary  Layer  Interaction  Investigations 
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the  constant  pressure  plateau  within  the  separated  region,  and  the  pressure  rise  to 
reattachment  as  described  in  reference  [19].  The  most  noticeable  aspect  of  the  pres¬ 
sure  profile  is  the  erratic  behavior  in  the  separated  region  when  the  Baldwin-Lomax 
transition  criteria  is  used.  Use  of  the  negative  shear  test  to  determine  the  transition 
location  results  in  a  much  smoother  pressure  profile  in  the  separated  region.  The  skin 
friction  profile  of  Figure  3.4  contains  several  regions  of  interest.  Using  the  Baldwin- 
Lomax  separation  criteria,  there  is  an  increase  in  the  friction  coefficient  leading  up 
to  the  sharp  drop  just  prior  to  separation.  The  Baldwin-Lomax  transition  criteria 
resulted  in  delayed  separation  while  the  negative  shear  test  criteria  led  to  a  closer 
match  with  the  experimental  data.  The  skin  friction  profile  beyond  reattachment 
shows  a  rapid  rise  to  the  ultimate  value  for  both  transition  criteria. 

The  solutions  obained  using  the  ATNSC  TVD  algorithm  clearly  demonstrate 
that  it  is  possible  to  obtain  very  accurate  estimations  of  separation  and  reattachment 
points  as  well  as  pressure  and  skin  friction  profiles.  ATNSC  solutions  are  obtained 
using  a  constant  CFL  number  of  0.95,  under  the  time  step  restriction  of  Eq  1.28.  The 
solution  is  monitored  until  no  change  is  observed  in  the  skin  friction  profile,  typically 
requiring  100000  time  steps  to  achieve  steady-state  convergence  with  =  2.2. 

3.4  High-Work  Low- Aspect-Ratio  Turbine 

The  ATNSC  scheme  is  next  apphed  to  a  rotor  cascade  whose  heat  transfer 
characteristics  were  studied  experimentally  by  Hippensteele,  Russell  and  Torres  [25] 
and  computationally  by  Boyle  [4].  The  comparisons  shown  herein  are  based  upon 
the  design  Reynold’s  number  of  7.6  x  10^  based  on  the  blade  chord. 

Figure  3.5  is  a  typical  C-type  grid  used  for  this  analysis.  The  actual  grid  used 
is  made  up  of  285  x  46  points.  Points  are  clustered  at  the  leading  and  traifing  edges 
for  improved  resolution.  25  points  are  placed  along  the  portion  of  the  C-type  grid 
representing  the  inlet. 

Consistent  with  subsonic  inflow  at  the  computational  inlet,  the  total  pressure 
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and  total  temperature  in  the  quiescent  region  infinitely  far  upstream  of  the  cascade 
are  required  as  boundary  conditions.  The  values  used  are  pt^  —  11.87  x  10'^ iV/m^ 
and  Tt^  =  303.00K.  The  static  pressure  at  station  3,  the  rotor  exit,  is  input  as 
the  exit  pressure.  In  particular,  pa  =  8.47  x  10^77/m^  For  the  viscous  calculations, 
the  surface  temperature  is  held  constant  at  =  333.337ir  consistent  with  Boyle’s 
work  [4,  5].  The  initial  conditions  applied  for  the  present  study  are  referred  to  as 
“cascade  tunnel  start”  conditions  because  of  the  analogy  to  the  starting  of  a  blow¬ 
down  cascade  tunnel.  The  domain  is  initialized  at  zero  velocity,  the  pressure  and 
temperature  corresponding  to  that  in  the  quiescent  region  upstream  of  the  inlet.  This 
is  analogous  to  placing  a  diaphragm  at  the  exit  of  the  computational  domain.  At 
time  to,  the  solution  is  started  and  a  centered  expansion  wave  propagates  upstream. 
It  is  also  possible  to  place  the  diaphragm  anywhere  in  the  computational  domain, 
but  placing  it  at  the  exit  avoids  the  formation  of  a  contact  surface  that  must  pass 
through  the  domain.  While  the  present  TVD  scheme  has  demonstrated  the  abihty 
to  resolve  such  a  contact  surface  in  very  fine  detail,  convergence  is  slowed  due  to 
the  fact  that  the  contact  surface  progresses  through  the  domain  at  the  convective 
velocity. 

When  the  cascade  tunnel  start  is  used  and  the  diaphragm  is  placed  at  the  exit  of 
the  computational  domain,  a  centered  expansion  wave  propagates  upsteam  through 
the  blade  passage  and  towards  the  inlet.  As  the  leading  edge  of  the  expansion  wave 
reaches  the  leading  edge  of  the  airfoil,  circulation  is  estabhshed  around  the  blade 
through  the  shedding  of  a  starting  vortex  from  the  airfoil.  The  vortex  is  convected 
downstream  and  eventually  exits  the  computational  domain  without  being  reflected. 

Figures  3.6  and  3.7  represent  the  solution  obtained  when  the  TVD  methodology 
is  apphed  to  this  case.  Figure  3.6  compares  the  surface  pressures  obtained  from  the 
Euler  TVD  solution  with  Boyle’s  Euler  solution  [4].  The  TVD  solution  shows  better 
agreement  with  the  data  except  in  the  area  of  the  trailing  edge.  This  is  explained 
by  the  fact  that  Boyle’s  calculations  were  performed  on  blade  with  a  cusp  added 
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to  the  trailing  edge.  The  current  solution  was  computed  for  a  blunt  trailing  edge, 
since  addinga  cusp  forces  the  location  of  the  stagnation  point  to  occur  at  the  sharp 
traihng  edge.  This  in  turn  specifies  the  circulation  around  the  blade. 

Figure  3.7  compares  the  Stanton  numbers  obtained  in  the  present  effort  with 
those  computed  by  Boyle  [4].  Boyle’s  turbulence  model  incorporated  the  effects 
of  pressure  gradient  and  freestream  turbulence  into  the  standard  Baldwin-Lomax 
model.  The  present  effort  utihzes  the  standard  Baldwin-Lomax  model  only.  The 
ATNSC  solution  follows  the  general  trend  of  the  experimental  data,  except  in  the 
region  just  downstream  of  the  transition  location  on  the  pressure  surface.  Boyle 
reports  that  his  calculations  predicted  early  transition  on  the  suction  surface.  It 
appears  that  his  calculations  actually  predict  suction  and  pressure  surface  transtion 
much  too  late.  The  ATNSC  calculations,  as  well  as  the  experimental  data  [25], 
show  that  the  large  increase  in  heat  transfer  near  the  suction  surface  traihng  edge 
is  actually  due  to  flow  separation,  not  transition.  This  is  shown  graphically  in 
Figure  3.8.  Transition  occurs  near  the  leading  edge  on  both  the  suction  and  pressure 
surfaces. 

The  steady-state  solution  obtained  with  the  TVD  formulation  compares  ex¬ 
tremely  well  against  the  experimental  data.  This  level  of  agreement  provides  an 
excellent  argument  for  the  use  of  TVD  schemes  in  computing  transonic  cascade 
flows. 

CFL  numbers  as  high  as  0.95  were  consistently  used  to  obtain  steady-state 
results.  In  fact,  the  CFL  number  was  dropped  to  0.5  only  if  a  contact  surface  was 
in  the  vicinity  of  the  rounded  traihng  edge  of  the  blade.  At  ah  other  times  the 
CFL  number  was  mcdntained  at  0.95.  This  is  in  contrast  to  CFL  numbers  as  low 
as  0.2  required  during  startup  and  only  as  high  as  0.8  to  maintain  stabihty  when 
using  the  MacCormack  scheme  [15,  17).  A  typical,  time-accurate  calculation  requires 
approximately  100,000  operator  sweeps  to  achieve  steady  state  convergence  for  the 
turbulent  computations. 
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Figure  3.9.  Flat  Plate  Mounted  in  Shock  Tube 

3.5  Unsteady  Shock-Induced  Heat  Transfer 

The  final  test  case  for  the  ATNSC  algorithm  is  the  prediction  of  unsteady  heat 
transfer  due  to  a  shock  wave  moving  down  a  flat  plate.  The  origin  of  this  test  case  is 
the  work  of  Smith  [39]  who  used  a  shock  tube  to  study  the  heat  transfer  to  a  sharp- 
edged  flat  plate,  creating  ratios  of  gas  temperature  to  surface  temperature  typical  of 
those  in  gas  turbine  engines.  A  schematic  of  Smith’s  experimental  appartus  is  shown 
in  Figure  3.9. 

The  case  under  consideration  here  is  representative  of  data  set  A  of  Smith  [39]. 
The  governing  parameters  are  the  shock  Mach  number  (Af,),  pressure  in  the  driven 
section  (Pi),  and  temperature  in  the  driven  section  (Ti).  The  wall  temperature 
on  the  fiat  plate  is  held  constant  in  the  calculations  at  Ti.  Using  M,  =  1.095, 
Pi  =  49102.800  NItvS,  and  Ti  =  297.428  K  results  in  a  shock  pressure  ratio  of  1.232, 
a  temperature  ratio  of  1.062,  a  steady  velocity  behind  the  shock  of  52.368  771/5,  a 
steady  Mach  number  behind  the  shock  of  0.147,  and  a  steady  Reynolds  number  be¬ 
hind  the  shock  of  92482.  Shock  Mach  number,  driven  section  pressure,  and  driven 
section  temperature  are  consistent  with  data  set  A  of  reference  [39].  Experimen¬ 
tal  measurements  of  the  shock  Mach  number  are  only  accurate  within  ±2%  and 
can  significantly  effect  the  level  of  agreement  [39]  between  theory,  experiment,  and 
numerical  solution.  Thus,  the  numerical  solutions  used  the  nominal  shock  Mach 
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Figure  3.10.  Grid  for  Heat  Flux  Solutions 


number  of  1.095  for  comparison.  A  solution  for  M,  =  1.117,  2%  above  nominal,  is 
also  presented. 

Initial  conditions  for  the  computations  consist  of  placing  the  shock  just  ahead  of 
the  plate  at  time  zero  by  establishing  pre-shock  and  post-shock  conditions  on  either 
side  of  the  point  selected.  The  shock  is  then  allowed  to  move  freely  as  time  progresses. 
With  the  temperature  held  at  Ti  the  normal  momentum  equation  is  solved  to  obtain 
the  pressure  at  the  wall,  in  combination  with  the  no-shp  constraint  at  the  waU.  The 
numerical  solution  is  sampled  at  a  point  5.080  x  10”^m  downstream  of  the  leading 
edge  to  obtain  the  heat  flux  at  this  point.  This  point  was  chosen  consistent  with 
the  first  sampling  point  of  Smith  in  the  shock  tube  experiment.  The  computations 
are  carried  out  to  a  time  of  approximately  one  miUisecond,  the  approximate  time 
of  transition  to  turbulent  flow  as  noted  by  Smith  [39] .  A  typical  grid  used  in  this 
study  is  shown  in  Figure  3.10.  This  grid  consists  of  201  points  in  the  axial  direction 
and  31  points  in  the  normal  direction.  The  grid  spacing  is  held  constant  in  the  axial 
direction  at  Aa:  =  2.540  X  10“^m  with  the  initial  spacing  from  the  wall  consistent 
with  =  2.2. 

Figure  3.11  is  the  solution  obtained  with  e  values  of  ei  =  £2  =  63  =  64  =  0.0, 
based  upon  the  numerical  experiments  of  the  author  [16].  Time  is  referenced  to 
the  time  of  the  shock  wave  passing  the  samphng  point.  The  peak  heat  flux  agrees 
well  with  experiment  for  both  shock  Mach  numbers.  The  Baldwin-Lomax  transition 
criteria  failed  to  initiate  transition  for  either  case.  Because  the  transition  mechanism 
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is  not  well  understood  for  this  experiment  [39],  the  boundary  layer  was  artificially 
forced  to  transition  at  the  leading  edge  in  the  computations.  At  the  higher  shock 
Mach  number  this  results  in  a  very  early  transition,  but  the  final  heat  flux  level  is  in 
good  agreement  with  the  thermocouple  data  after  approximately  1.8  msec.  At  the 
lower  shock  Mach  number,  transition  again  occurs  early  but  the  final  heat  flux  level 
agrees  well  with  the  heat  flux  gage  data.  It  is  clear  from  these  results  that  much 
more  effort  is  needed  in  an  attempt  to  understand  the  transition  mechanism  for  this 
situation.  AU  ATNSC  solutions  are  undertaken  at  a  CFL  of  0.95,  with  the  time  step 
restriction  of  Eq  1.28.  This  results  in  approximately  25,000  sweeps  to  arrive  at  a 
time  of  2.0  msec. 

3.6  Conclusions 

TVD  methodolgy  has  been  applied  to  problems  not  previously  examined  using 
TVD  schemes.  Early  TVD  research  concentrated  mainly  on  supersonic  and  hyper- 
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Figure  3.11.  Heat  Flux  History 
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sonic  flows,  both  inviscid  and  viscous,  and  was  almost  solely  directed  toward  obtain¬ 
ing  steady-state  solutions.  Later  effort  by  the  author  [16]  extended  the  TVD  method¬ 
ology  to  inviscid  transonic  cascade  flows,  viscous  flows  with  shock-induced  laminar 
boundary-layer  separation,  and  unsteady  laminar  flows  with  signifcant  shock-induced 
heat  transfer.  The  current  effort  has  extended  the  TVD  methodology  to  both  steady 
and  unsteady  turbulent  flows. 

Shock-boundary-layer  interaction  has  been  studied  by  numerous  researchers 
using  highly  regarded  algorithms  such  as  the  MacCormack  [29],  Dawes  [10],  Beam- 
Warming  [11],  and  Newton  [28]  methods.  While  acceptable  predictions  of  the  pres¬ 
sure  profile  in  the  boundary-layer  have  often  been  computed,  TVD  schemes  have 
enabled  researchers  to  accurately  compute  the  skin  friction  profile  [16].  The  ATNSC 
algorithms  finally  provide  the  means  to  accurately  compute  pressure  profiles,  sepa¬ 
ration  and  reattachment  locations,  and  skin  friction  profiles  in  good  agreement  with 
the  available  experimental  data.  Results  given  in  Section  3.3  are  testament  to  this. 

Cascade  flows  are  currently  of  great  interest  to  gas  turbine  engine  designers 
and  researchers  [2,  14,  13].  Analysis  of  these  flows  is  a  severe  test  of  an  algorithm 
because  of  the  wide  Mach  number  range,  typically  0.3  <  M  <  1.3  ,  and  the  fact 
that  the  flow  is  confined  in  a  passage  where  wave  systems  tend  to  reflect  back  into 
the  domain.  Results  presented  in  Section  3.4  show  that  the  TVD  schemes  of  Sec¬ 
tions  2.2.1  and  2.2.2  yield  steady-state  results  in  good  agreement  with  experiment, 
compared  to  other  methods. 

Unsteady  shock-induced  heat  transfer  has  been  studied  theoretically  [30,  31], 
experimentally  [14,  39],  and  has  recently  become  of  interest  computationally  [16]. 
The  increased  interest  is  due  to  the  enhanced  system  performance  available  through 
accurate  knowlege  of  the  heat  transfer  [6].  However,  it  is  not  uncommon  for  computed 
heat  flux  values  to  be  an  order-of-magnitude  different  from  experimental  values  [18]. 
The  ATNSC  algorithms  represent  a  significant  advancement  in  the  state-of-the-art 
for  computing  shock-induced  heat  transfer. 
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5.7  Further  Research 


Performance  of  the  ATNSC  ailgorithms  with  enhanced  turbulence  models,  ac¬ 
counting  for  freestream  turbulence,  pressure  gradient,  etc.,  should  be  investigated. 
These  enhancements  were  not  possible  in  this  hmited  first  evaluation  of  TVD  per¬ 
formance  with  turbulence  modeling  included.  The  current  investigation  shows  that 
implementation  of  the  simple  Baldwin-Lomax  algebraic  turbulence  model  correctly 
predicts  both  skin  friction  and  heat  flux  in  the  turbulent  region. 

Overall,  the  TVD  based  viscous  algorithms  perform  exceptionally  well  on  the 
test  cases  herein.  Emphasis  must  be  placed  on  applying  these  algorithms  to  even 
more  rigorous  test  cases  so  as  to  gcdn  an  even  greater  understanding  of  their  weak¬ 
nesses  as  well  as  their  strengths. 
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